function [A, b] = process_boundary_2d(A, b, mesh, pde)
%% Robin
matrix_A = sparse(mesh.Nb,mesh.Nb);
vector_b = zeros(mesh.Nb,1);
cr = @(x,y) pde.cor(x,y).*pde.r(x,y);
for k = 1:size(pde.bdary,2)
    if pde.bdary(1,k) == 2
        bde = boundaryedges(mesh, k);
        for n = 1:size(bde,1)
            points = mesh.node(bde(n,2:3),:);
            AB = [points(2,1)-points(1,1); points(2,2)-points(1,2)];
            vector_n = [0, 1; -1, 0] * AB;
            vector_n = vector_n./norm(vector_n);
            cq = @(x,y) pde.cor(x,y).*(vector_n(1)*pde.dudx(x,y) + vector_n(2)*pde.dudy(x,y) + pde.r(x,y).*pde.u(x,y));
            nk = bde(n,1);
            E = mesh.node(mesh.elem(nk,:),:);
            V = E(1:3,:);
            [gauss_points_loc, gauss_weights_loc] = gauss_lengdre_loc(points, "1d");
            psi = basis_function_loc(V, gauss_points_loc(:,1), gauss_points_loc(:,2), mesh.basis_type);
            fcq = cq(gauss_points_loc(:,1), gauss_points_loc(:,2));
            fcr = cr(gauss_points_loc(:,1), gauss_points_loc(:,2));
            for i = 1:mesh.Nlb
                f = fcq .* psi{i};
                r = gauss_weights_loc * f;
                vector_b(mesh.elem(nk,i),1) = vector_b(mesh.elem(nk,i),1) + r;
                for j = 1:mesh.Nlb
                    f = fcr .* psi{i} .* psi{j};
                    r = gauss_weights_loc * f;
                    matrix_A(mesh.elem(nk,i),mesh.elem(nk,j)) = matrix_A(mesh.elem(nk,i),mesh.elem(nk,j)) + r;
                end
            end
        end
    end
end
A = A + matrix_A;
b = b + vector_b;

%% Neumann
vector_b = zeros(mesh.Nb,1);
for k = 1:size(pde.bdary,2)
    if 1 == pde.bdary(1,k)
        bde = boundaryedges(mesh, k);
        for n = 1:size(bde,1)
            points = mesh.node(bde(n,2:3),:);
            AB = [points(2,1)-points(1,1); points(2,2)-points(1,2)];
            vector_n = [0, 1; -1, 0] * AB;
            vector_n = vector_n./norm(vector_n);
            cp = @(x,y) pde.cor(x,y).*(vector_n(1)*pde.dudx(x,y) + vector_n(2)*pde.dudy(x,y));
            nk = bde(n,1);
            E = mesh.node(mesh.elem(nk,:),:);
            V = E(1:3,:);
            [gauss_points_loc, gauss_weights_loc] = gauss_lengdre_loc(points, "1d");
            psi = basis_function_loc(V, gauss_points_loc(:,1), gauss_points_loc(:,2), mesh.basis_type);
            fcp = cp(gauss_points_loc(:,1), gauss_points_loc(:,2));
            for i = 1:mesh.Nlb
                f = fcp .* psi{i};
                r = gauss_weights_loc * f;
                vector_b(mesh.elem(nk,i),1) = vector_b(mesh.elem(nk,i),1) + r;
            end
        end
    end
end
b = b + vector_b;

%% Dirichlet
if ~isempty(mesh.mesh_type)
    bdn = mesh.bdn;
    A(bdn,:) = 0;
    A(sub2ind(size(A), bdn, bdn)) = 1;
    b(bdn) = pde.u(mesh.node(bdn,1),mesh.node(bdn,2));
else
    for k = 1:size(pde.bdary,2)
        if pde.bdary(1,k) == 0
            bdn = boundarynodes(mesh, k);
            A(bdn,:) = 0;
            A(sub2ind(size(A), bdn, bdn)) = 1;
            if pde.err
                b(bdn) = pde.u(mesh.node(bdn,1),mesh.node(bdn,2));
            else
                b(bdn) = pde.D(mesh.node(bdn,1),mesh.node(bdn,2),k);
            end
        end
    end
end
end

%% boundarynodes
function bdn = boundarynodes(mesh, index)
switch index
    case 1
        if "P2" == mesh.basis_type
            bdn = mesh.bdn(1:2*mesh.Nx+1);
        else
            bdn = mesh.bdn(1:mesh.Nx+1);
        end
    case 2
        if "P2" == mesh.basis_type
            bdn = mesh.bdn(2*mesh.Nx+1:2*mesh.Nx+2*mesh.Ny+1);
        else
            bdn = mesh.bdn(mesh.Nx+1:mesh.Nx+1+mesh.Ny);
        end
    case 3
        if "P2" == mesh.basis_type
            bdn = mesh.bdn(2*mesh.Nx+2*mesh.Ny+1:4*mesh.Nx+2*mesh.Ny+1);
        else
            bdn = mesh.bdn(mesh.Nx+1+mesh.Ny:2*mesh.Nx+1+mesh.Ny);
        end
    case 4
        if "P2" == mesh.basis_type
            bdn = mesh.bdn(4*mesh.Nx+2*mesh.Ny+1:end);
        else
            bdn = mesh.bdn(2*mesh.Nx+1+mesh.Ny:end);
        end
        bdn = [bdn; 1];
    otherwise
        error("Invalid edge number");
end
end

%% boundaryedges
function bde = boundaryedges(mesh, index)
switch index
    case 1
        bde = mesh.bde(1:mesh.Nx, 1:3);
    case 2
        bde = mesh.bde(mesh.Nx+1:mesh.Nx+mesh.Ny, 1:3);
    case 3
        bde = mesh.bde(mesh.Nx+mesh.Ny+1:2*mesh.Nx+mesh.Ny, 1:3);
    case 4
        bde = mesh.bde(2*mesh.Nx+mesh.Ny+1:end, 1:3);
    otherwise
        error("Invalid edge number");
end
end